Superfluid to solid crossover in a rotating Bose-Einstein condensed gas 
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The properties of a rotating Bose-Einstein condensate confined in a prolate cylindrically symmetric 
trap are explored both analytically and numerically. As the rotation frequency increases, an ever 
greater number of vortices are energetically favored. Though the cloud anisotropy and moment 
of inertia approach those of a classical fluid at high frequencies, the observed vortex density is 
consistently lower than the solid-body estimate. Furthermore, the vortices are found to arrange 
themselves in highly regular triangular arrays, with little distortion even near the condensate surface. 
These results are shown to be a direct consequence of the inhomogeneous confining potential. 
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One of the most striking properties of liquid '*He(II) 
is its ability to mimic the behavior of a solid body when 
subjected to uniform rotation. Since the superfluid ve- 
locity field Vs is irrotational (V x Vg = 0), the superfluid 
component of '*He(II) might be expected to remain at 
rest while the normal component rotates with the con- 
tainer. In fact, for sufficiently large values of the rotation 
frequency Q, the entire fluid is found to rotate like a clas- 
sical liquid at all temperatures [Q . The paradox may be 
resolved by assuming that the superfluid is threaded by 
quantized vortices. These are singularities of v^, around 
which the phase of the superfluid order parameter in- 
creases by 27r. Although the mechanisms for the spin-up 
of the superfluid are not fully understood, at equilibrium 
the vortices must flow with the normal velocity due to 
the the mutual friction between superfluid and normal 
components In addition to considerable indirect evi- 
dence for this hypothesis, small numbers of vortices have 
been imaged directly in rotating superfluid ^He 

In this work we show how the presence of quantized 
vortices can allow a Bose-Einstein condensate (BEC) to 
mimic a classical fluid under rotation, as has been sug- 
gested by recent experiments at JILA In these ex- 
periments, a trapped gas of ultracold ^^Rb atoms is spun 
up, and then cooled through the Bose-Einstein conden- 
sation transition. For small values of f2, the condensate 
density is found to assume its usual non-rotating shape, 
while the thermal cloud bulges outward. This corrobo- 
rates previous evidence that the condensate behaves as 
an irrotational superfluid The condensate density 

profile undergoes a sudden change at a value of Q that is 
comparable to the thermodynamic critical frequency for 
the stability of a single vortex As fl increases fur- 
ther, the shape of the condensate gradually approaches 
that of the thermal cloud. This suggests that for any 
given value of and temperature, the condensate con- 
tains the appropriate number and distribution of vortices 
for thermodynamic equilibrium. In contrast, when no 
appreciable thermal fraction is present, higher rotation 
frequencies are generally required to nucleate vortices in 



BECs 

We present two key results, which also bear on issues 
raised by recent experiments at MIT First, the vor- 
tices are arranged in extremely regular triangular arrays: 
even near the condensate surface, little circular distor- 
tion is found. Second, the number of vortices is con- 
sistently lower than that required to ensure solid-body 
rotation throughout the condensate. 

To make explicit comparison with the recent JILA ex- 
periment, we consider the case of iV = 200, 000 atoms 
of ^^Rb, confined in a cylindrically symmetric trap with 
radial frequency uip /27r = 8 Hz, and anisotropy A = 
LUz/ujp = |. Unless stated explicitly, our units of energy, 
angular frequency, length, and time are given by hujp, 

respectively, 



dp = y^h/AIujp 3.845 /im, and lu 
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where M is the atomic mass and h is Planck's constant h 
divided by 27r. We work in a frame that rotates with an- 
gular frequency ft about the z axis. The time-dependent 
Gross-Pitaevskii (GP) equation which governs the 
dynamics of the condensate wavefunction tp oi a, dilute 
BEC at zero temperature, is then given by: 



iftV(r, t)^[T + Ftrap + Vk- ni,] iPir, t), 



(1) 



with kinetic energy T = — iV^, trap potential Vtrap — 
i (p^ -I- A^z^), and angular momentum component = 
i{ydx — xdy). The effects of atomic interactions are 
included in the nonlinear term Vh = 47r?7|V'P, rj = 
Na/dp, where a — 5.29 nm is the scattering length for 
®^Rb collisions [Q. We use the normalization condition 
J dr\'ip{r,t)\'^ — 1. In equilibrium in the rotating frame, 
ip{v,t) = e~'^*7/'(r), where n is the chemical potential. 

To estimate the properties of a rotating condensate, 
such as the aspect ratio and the number of vortices, we 
consider two tractable cases: a single vortex applicable 
for small SI, and multiple vortices relevant to high f2 
where the condensate is expected to behave essentially 
as a rigid body. 

With one vortex at the center of the trap, i/i = 
l^/ije*"^, where (p is the polar angle. In the large-TV 
or Thomas- Fermi (TF) limit, the condensate density is 
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n ) / 4:717] when that quan- 



tity is positive, and is zero elsewhere. [Q. The inner 
cutoff defines the vortex core size or the healing length 
^; for /X 3> 1, one obtains ^ w l/y/2jl ^ l/Ro, where 
i?o = {15r]Xy^^ is the TF radius along p in the ab- 
sence of a vortex. A straightforward calculation shows 
that, for large Rq, the TF radius for an isolated vortex 
is Rp « Ro[l + (3/2i?;])ln(2i?o/0] and the condensate 
aspect ratio is AJpp = Rp/Rz ~ A[l + (l/2i?§)]. Assuming 
this result depends only weakly on vortex position, and is 
additive with respect to the number of vortices Ny , then 
in explicit units A^pp ~ A[l + ^Ny{dp/ Ro)"^] at larger ft. 

For large N^, as we will show below, the condensate 
rotates almost as a solid body, so the rotating-frame ve- 
locity operator — — iV — fii x r can be neglected. 
1 

2 "r 2 

ens the radial potential, Vtrap ^ ^(1 

r!2)3/io and A^j, 



Since T — flL^ — — ^il^p^, rotation effectively soft- 



In 



this case, Rp = Ro/{l - ^2)3/10 ^nd A^^ = A/Vl - f^^. 
The number of vortices is the line integral of the phase 
gradient around the cloud perimeter; assuming the solid- 
body value of the tangential velocity ^Rp, then the areal 
vortex density is — Nv/iiR^p — JI/tt, and Nf" — 

If the vortices form a regular array at large f2, then 
the lattice constant h should be comparable to the aver- 
age separation between vortices ~ ^Ji: jVl. For a 
triangular array centered at the origin |l^ , p^ , the vor- 
tices arrange themselves in concentric hexagonal rings 
labelled by ring index r, such that — I + 3r(r -I- 1). 
Assuming the superfluid velocity exactly matches the 
solid-body value midway between nearest-neighbor vor- 
tices (where the two nearest rotational fields exactly can- 
cel), then iV„ = ri&2(r -I- \Y and h « y'S/fl for large r. 
Since Jlmax = 1 in a harmonic trap, the smallest vor- 
tex separation is &i„in ~ V^dp in explicit units. When 
the vortex cores begin to overlap significantly (6 ~ C)j 
the system might undergo a phase transition, possibly 
into a state akin to a quantum Hall insulator |15| ; since 
^(p = 0, 2 = 0) = = l/^o(l - ^"^f"^ in the TF 
limit, the value of f2 for this to occur must become ex- 
tremely close to unity: 1 — Rq^ ■ 

The stationary solutions of the GP equation in the 
rotating frame, defined as local minima of the free en- 
ergy {E) = fiN — ^{Vh), are found numerically by norm- 
preserving imaginary time propagation using an adap- 
tive stepsize Runge-Kutta integrator. The wavefunc- 
tion is solved on a three-dimensional Cartesian mesh 
within a discrete-variable representation [p^ based on 
Gauss-Hermite quadrature, and is assumed to be even 
under reflection in the z = plane. The initial state is 
taken to be the TF wavefunction with a phase <I>(x, y) — 
J2xo,yo tan~M(2/ - yo)/ix - xq)], where {xq, yo) are vor- 
tex positions in a regular array centered at the origin. 
The GP equation for a given value of is propagated in 
imaginary time until the fluctuations in both p and the 



norm become smaller than 10~^^. The condensate densi- 
ties integrated down x and z are then flt to a TF profile 
using a nonlinear least-squares analysis, where densities 
lower than 0.1% of the maximum value are discarded. For 
the vortex-free condensate, this yields an aspect ratio of 
0.645, which is 3% larger than the TF value of |. 

The resulting equilibrium configurations are sensitive 
to the initial vortex distributions. Fig. |l| shows three 
different solutions of the GP equation (|^) for Q = 0.45. 
These were obtained using seed arrays with rhombohe- 
dral (left), square (center), and triangular (right) sym- 
metries, respectively. Though observables such as the 
energy, angular momentum, and cloud aspect ratio are 
all comparable, they each have different vortex numbers 
and arrangements. Though a complete survey of possible 
configurations is beyond the scope of the present work, 
for all cases considered the initial rhombohedral vortex 
distribution is found to yield the final state with lowest 
energy; for larger Q this symmetry gives rise to equilib- 
rium arrays that are generally triangular (see below). 
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FIG. 1. Three stationary states are shown for f2 = 0.45. 
From left to right, the values {E /N , {L^) /N} in 
scaled units are {7.20,11.28,2.31}, {7.21,11.29,2.33}, and 
{7.24, 11.28, 2.19}, respectively. 

The central results of the present work are shown in 
Fig. ^, which depicts equilibrium solutions for 0.25 < 
f2 < 0.95. A single vortex at the origin has appeared 
hy — 0.35; the thermodynamic critical frequency (the 
energy difference between states with zero and one vor- 
tex, divided by Ti) is fic = 0.30. This value is slightly 
lower than the experimental value 0.32 < Q,c < 0.38; 
since Vtc ^ N^'^^^, perhaps there are fewer atoms in the 
condensate at vortex nucleation. (The dynamic critical 
frequency, at which the first collective mode becomes neg- 
ative, is somewhat higher: fli, — 0.46). With a vortex, 
the cloud aspect ratio changes to A = 0.663; using the 
fitted values for the nonrotating cloud A = 0.645 and 
Ro = 4.86, the TF prediction is A^p = 0.659. As Q. 
continues to increase, so does the aspect ratio; the cloud 
becomes spherical for 0.75 < < 0.8 (consistent with 
the experimental results) and highly oblate for fl = 0.95, 
at which A' = 1.8. 

As shown in Fig. ||, the solid-body estimate of the cloud 
anisotropy Ag^^ tracks (but consistently exceeds) our nu- 
merical values; in contrast, A^^p is always too small. For 
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example, when J7 = 0.95 one obtains A^^^ = 2.00 and 
A^TF = 1.50 with Ny = 65 (see Fig. |). The number of 
enclosed vortices is not known a priori, however; using 
the sohd-body estimate Nf" = 89 for n = 0.95 yields 
the much improved Al^p = 1.83. Another indication 
that the condensate is behaving classically at large is 
the moment of inertia / (inset of Fig. ^). The effective 
value / — (Lz)/^ is always lower than the solid-body 
/ = (^2 + 2/2)^ but is within 4% hy n ^ 0.95. 
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FIG. 2. Condensate densities integrated along x (upper 
row) and z (lower row) are shown for Q — 0.25, 0.35, 0.55, 
and 0.65 (first data set, left to right) and Cl — 0.75, 0.80, 
0.90, and 0.95 (second data set, left to right). Each frame is 
20dp ~ 77 on a side. 



The number of vortices at equilibrium is always con- 
siderably lower than the solid-body prediction, as in pre- 
vious experimental observations Since the numeri- 
cal solutions are stationary in the rotating frame, this 
discrepancy cannot be explained by positing that the 
vortex array rotates more slowly than the trap. Con- 
sider the cases ft — 0.55, 0.8, 0.9, and 0.95 shown in 
Fig. ^ which approximate centered triangular arrays 
with Ny — 1 + 3r(r -|- 1), r = 1 — 4, respectively. The 
average vortex spacing is found to follow the prediction 
b = y/S/fl to within 3%. An additional hexagonal ring of 
vortices could therefore fit comfortably within the cloud. 
For = 0.95, 5b = 8.89 is smaller than the radius 
Rp — 9.41, and r = 5 corresponds to = 91 which 
is close to the solid-body prediction N!^^ — 89. For 
Ny = 169 (r = 7), which is comparable to the largest 
array in experiments at MIT ||^ , the missing = 8 ring 
implies that the equilibrium number of vortices is of order 



20% lower than the solid-body prediction. 

The absence of the last ring might be due to the fact 
that vortices in this low-density region would significantly 
overlap because of their large core size. Assuming that 
the vortex diameter is twice the local healing length, then 

with ^{p,z = 0) = l/Rp^{l-n^){l- p'^/Rl) one ob- 
tains a critical vortex displacement pc ^ 9 for Q = 0.95. 
In fact, the energy of a uniform array of vortices in a ro- 
tating cylinder is also minimized if there exists a 'vortex- 
free strip' the size of approximately one ring near the edge 
of the vessel @, i.e. Ny = 2nR'^n/K-S, where S ~ N''^. 
This correction is due to the contribution to the energy 
of strictly irrotational flow in the region between the last 
vortex and the superfluid surface. 
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FIG. 3. The cloud aspect ratio A = Rp/Rz and moment 
of inertia / (inset) are shown as a function of rotation fre- 
quency Q. The numerical results (solid lines) for A are ob- 
tained by a TF fit to the cloud profile, and / = {Lz)/NQ.. 
The solid-body (dashed lines) result for the cloud anisotropy 
corresponds to — Q?, while the corresponding moment 

of inertia is given by / = {x^ + y^). 

The existence of a vortex-free region in trapped con- 
densates is confirmed by evaluating the change in con- 
densate phase around a contour in the a;y-plane of in- 
creasing radius R from the origin. This is accomplished 
by calculating the spatial derivatives of the numerical 
data in order to determine v = V$, interpolating the 
results onto a one-dimensional azimuthal grid with 2000 
points, and evaluating the line integral §v ■ d\ numeri- 
cally using a trapezoidal rule. The results for 11 = 0.75 
and 0.95 are given in Fig. I On average, the number of 
vortices follows the solid-body expression U,R^ for small 
rings, but begins to lag noticeably as R ^ Rp even be- 
fore the vortex-free strip is reached. The velocity field 
for the = 0.95 case, shown in Fig. ||, is small in the 
rotating frame everywhere except for the rotational cur- 
rents near the vortex cores and the irrotational flow near 
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the surface. 
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FIG. 4. The number of vortices within a circular contour 
centered at the origin are shown as a function of radius R 
(solid lines) for Q, — 0.75 and 0.95. The solid-body predictions 
Q,R^ (dashed lines) are shown for comparison. The vertical 
dotted lines denote the TF fit for the radial radius. 

In order to further explore this issue, consider a model 
wavefunction with constant amplitude and phase given 

by ^ix,y) = J2xo.yQ^^'^~^iiy - yo)/{x - ^o)], where 
ixo,yo) are vortex positions in a centered triangular ar- 
ray with lattice constant b. For Ny = 61 (r = 4), 
the vortex velocities v — |V<i>| on successive hexago- 
nal rings rir are v — -^{3.63, 7.23, 10.69, 13.57}. Since 
v{nr = 4) < 4:v{nr = 1) by 7%, the angular velocity of 
the last ring cannot attain the solid-body value for any 
choice of b. For large arrays, this mismatch in veloci- 
ties varies as {R/Rpf , and is why significant distortion 
of the vortex array from triangular is expected near the 
superfluid surface [O. 




FIG. 5. The velocity field v in the a;j/-plane is represented 
by arrows for the = 0.95 case. The left and right images 
correspond to the lab and rotating frames, respectively. 

The question that immediately arises is: why are the 
vortex arrays observed in confined condensates so per- 
fectly triangular, even very near the surface? One pos- 
sible explanation is that a displaced vortex will precess 
around the origin even in the absence of other vortices, 
due to the inhomogeneous external potential. Neglecting 
vortex curvature (which from Fig. ^ is evidently negli- 
gible at large il), the additional contribution to the ve- 



locity is ^; = [R/{Rl - R^)] ln(C/i?p) in the TF limit 0. 
Let us return to the case considered above, with r — 
and choose fl — 0.95 for concreteness. Assuming Rp — 
Ro/{l - 1^2)3/10 and imposing 3.63/6 + v{R = b) = fib, 
one obtains b = 1.98 and v ^ {1.88, 3.76, 5.62, 7.37}. 
Thus, including the effect of precession, the solid-body 
value i; = 4xl.88ati? = 6 now exceeds the velocity of 
the last ring 7? = 46 by only 2%. 

In conclusion, we have explored the crossover of a con- 
fined Bose-Einstein condensate from that of an irrota- 
tional superfluid to a solid body with increasing rotation. 
The external potential is shown to strongly influence the 
density and arrangement of the resulting vortices. Many 
related issues remain unresolved, however, among them 
the spin-up of the superfluid by the thermal cloud, the 
upper critical frequency, and the approach to a quantum 
Hall state; these will be the subject of future work. 
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